
C     A COLLECTION OF COMMON ROUTINES, MOSTLY FOR SETUP OF THE
C     TRAVELLING SALESMAN PROBLEM.
C     DOUBLE PRECISION VERSION
C
C  SETPTS:  COMPUTE COORDINATES OF VERTICES
C
C  THIS ROUTINE TAKES AS ARGUMENTS THE PARAMETERS OF THE CRYSTAL
C  AND THE DESIRED RANGES OF MILLER INDICES AND COMPUTES THE COORDINATES
C  OF THE DETECTOR CORRESPONDING TO THE DESIRED REFLECTIONS.
C
C  INPUT PARAMETERS:
C     ORIENT  MATRIX OF VECTORS DEFINING RECIPROCAL LATTICE
C     LAMBDA  WAVELENGTH OF X-RAY BEAM IN INVERSE ANGSTROMS
C     HLO, HHI, KLO, KHI, LLO, LHI  RANGES OF MILLER INDICES H, K, L
C
C  OUTPUT:
C     NO PARAMETERS ARE CHANGED.  IN COMMON /POINTS/, N IS SET TO THE
C     NUMBER OF POINTS TO WHICH IT IS POSSIBLE TO MOVE THE DETECTOR,
C     AND THE COORDINATES OF THOSE POINTS ARE PLACED IN PHI, CHI,
C     AND TWOTH.
C
      SUBROUTINE SETPTS (ORIENT, LAMBDA, HLO, HHI, KLO, KHI, LLO, LHI)
C
      COMMON /MSTPRM/ MSGLVL
      COMMON /POINTS/ PHI, CHI, TWOTH, N
      INTEGER N, MSGLVL
      INTEGER INC ,INC2, KSTART,LSTART,KEND,LEND,MOD
      DOUBLE PRECISION    PHI(20000), CHI(20000), TWOTH(20000)
C
      DOUBLE PRECISION ORIENT(3,3), LAMBDA
      INTEGER HLO, HHI, KLO, KHI, LLO, LHI
      INTEGER H, K, L, IMPOSS
      DOUBLE PRECISION P, C, T
      LOGICAL POSIBL
      DOUBLE PRECISION DMIN1,DMAX1
      DOUBLE PRECISION MINP,MINC,MINT,MAXP,MAXC,MAXT
      MINP =  360.0D+0
      MINC =  180.0D+0
      MINT =  155.0D+0
      MAXP =    0.0D+0
      MAXC = -180.0D+0
      MAXT =  -55.0D+0
C
      IF (MSGLVL .GE. 2) WRITE (*,*) 'ENTERING SETPTS'
      N=0
      IMPOSS = 0
      DO H = HLO, HHI
        IF (MOD(H-HLO,2).EQ.0) THEN
          INC = 1
          KSTART = KLO
          KEND = KHI
          ELSE
            INC = -1
            KSTART = KHI
            KEND = KLO
          ENDIF
        DO K = KSTART, KEND, INC
          IF (MOD((H + K - HLO -KLO),2).EQ.0) THEN
            INC2 = 1
            LSTART = LLO
            LEND = LHI
            ELSE
              INC2 = -1
              LSTART = LHI
              LEND =LLO
            ENDIF
          DO L = LSTART, LEND, INC2
            IF (MSGLVL .GE. 2) WRITE (*,*) 'POINT H, K, L =', H, K, L
            IF (MSGLVL .GE. 2) WRITE (*,*) 'ENTERING ANGLES'
            IF (MSGLVL .GE. 2) WRITE (*,*) 'LAMBDA =', LAMBDA
            CALL ANGLES (H, K, L, ORIENT, LAMBDA, 0.0D0, P, C, T, 
     $          POSIBL)
            IF (POSIBL) GOTO 90
                IMPOSS = IMPOSS + 1
                IF (MSGLVL .GE. 2) WRITE (*,*) 'POINT H, K, L =',
     $              H, K, L, ' IS NOT POSSIBLE'
                GOTO 100
90            CONTINUE
                N = N+1
                PHI(N) = P
                CHI(N) = C
                TWOTH(N) = T
                MINP = DMIN1(MINP,P)
                MINC = DMIN1(MINC,C)
                MINT = DMIN1(MINT,T)
                MAXP = DMAX1(MAXP,P)
                MAXC = DMAX1(MAXC,C)
                MAXT = DMAX1(MAXT,T)
                IF (MSGLVL .GE. 2) WRITE (*,*) N, 'H,K,L=',
     $              H, K, L, 'PHI,CHI,TWOTH=', P, C, T
c********
                WRITE (*,*) n, p, c, t
c********
          END DO
        END DO
      END DO
100   CONTINUE
      IF (MSGLVL .GE. 1)
     $    WRITE(*,*) IMPOSS, ' POINTS WERE NOT POSSIBLE.'
      IF (MSGLVL .GE. 1) WRITE (*,*) N, 'POINTS TO VISIT.'
      WRITE(*,*) 'PHI RANGE:   ',MINP,MAXP
      WRITE(*,*) 'CHI RANGE:   ',MINC,MAXC
      WRITE(*,*) 'TWOTH RANGE: ',MINT,MAXT
      IF (MSGLVL .GE. 2) WRITE (*,*) 'LEAVING SETPTS'
      RETURN
      END
C  TCOST:  RETURN THE TOTAL COST OF A TOUR
C
C  INPUT:
C      TOUR    TOUR(V) IS THE VERTEX FOLLOWING V IN THE TOUR.
C              THE POINTS ARE IN COMMON /POINTS/ AS DESCRIBED IN SETPTS.
C
      DOUBLE PRECISION FUNCTION TCOST (TOUR)
C
      COMMON /MSTPRM/ MSGLVL
      COMMON /POINTS/ PHI, CHI, TWOTH, N
      INTEGER N, MSGLVL
      DOUBLE PRECISION    PHI(20000), CHI(20000), TWOTH(20000)
C
      EXTERNAL COST
      DOUBLE PRECISION COST
      INTEGER TOUR(N)
      INTEGER V, W, I, HOURS, MINITS, SECNDS
      DOUBLE PRECISION C
C
      IF (MSGLVL .GE. 2) WRITE (*,*) 'ENTERING TCOST'
      V = 1
      TCOST = 0.0D0
      DO 100 I = 1, N
          W = TOUR(V)
          C = COST(V,W)
          IF (MSGLVL .GE. 3)
     $        WRITE (*,*) 'EDGE:', V, W, 'COST:', C
          TCOST = TCOST + C
          V = W
          IF (V .EQ. 1  .AND.  I .NE. N) GOTO 800
100       CONTINUE
      IF (V .NE. 1) GOTO 800
      HOURS = TCOST / 3600
      MINITS = (TCOST-3600*HOURS) / 60
      SECNDS = TCOST - 3600*HOURS - 60*MINITS
      IF (MSGLVL .GE. 1)
     $    WRITE (*,*) 'COST OF TOUR IS', TCOST, ' SEC.'
      IF (MSGLVL .GE. 1)
     $    WRITE (*,*) HOURS, 'HR.', MINITS, 'MIN.', SECNDS, 'SEC.'
      IF (MSGLVL .GE. 2) WRITE (*,*) 'LEAVING TCOST'
      RETURN
C
800   IF (MSGLVL .GE. 1)
     $    WRITE (*,*) 'TCOST DETECTS ILLEGAL TOUR'
      STOP
      END
C  ANGLES:  GIVEN MILLER INDICES OF A REFLECTION,
C           COMPUTE POSITIONING INFORMATION FOR
C           THE DETECTOR.
C
C  FROM MATT SMALL, APRIL 5, 1984.
C
C  INPUT PARAMETERS:
C      IH, K, L    MILLER INDICES
C      ORIENT      3 BY 3 MATRIX OF VECTORS DEFINING THE RECIPROCAL
C                  LATTICE
C      LAMBDA      WAVELENGTH OF X-RAY BEAM (IN INVERSE ANGSTROMS)
C      OMEGA       'MUST BE KEPT AT 0.0' - FINN NIELSEN
C
C  OUTPUT PARAMETERS:
C      FI, KHI, TWOT  CALCULATED ANGLES PHI, CHI, AND TWO*THETA
C      POSIBL         FALSE IF IT IS IMPOSSIBLE TO MOVE
C                     TO THE REFLECTION.
C
C
      SUBROUTINE ANGLES (IH, K, L, ORIENT, LAMBDA, OM,
     $                   FI, KHI, TWOT, POSIBL)
      INTEGER IH, K, L
      DOUBLE PRECISION    FI, KHI, TWOT, LAMBDA, OM, OMEGA, ORIENT(3,3)
      LOGICAL POSIBL
C
      DOUBLE PRECISION PI, RH, RK, RL, COSOMG, X, Y, Z, D, DUM
      DOUBLE PRECISION T1, T2, T3, SINKHI, Q, R, COSFI, SINFI
      DATA PI /3.14159265368979/
C
      POSIBL = .TRUE.
      OMEGA = OM/180.0*PI
      RH = IH
      RK = K
      RL = L
      COSOMG = DCOS (OMEGA)
      X = ORIENT(1,1)*RH + ORIENT(1,2)*RK + ORIENT(1,3)*RL
      Y = ORIENT(2,1)*RH + ORIENT(2,2)*RK + ORIENT(2,3)*RL
      Z = ORIENT(3,1)*RH + ORIENT(3,2)*RK + ORIENT(3,3)*RL
      D = DSQRT (X*X + Y*Y + Z*Z)
      DUM = LAMBDA*D/2.
      IF (DUM .LT. .000001 .OR. DUM .GE. 1.0) GOTO 7
      TWOT = DASIN(DUM)*2.
      T1 = X/D
      T2 = Y/D
      T3 = Z/D
      IF (DABS(T3) .LT. DABS(COSOMG)) GOTO 2
7         POSIBL = .FALSE.
          GOTO 998
2     SINKHI = 0.0 - T3/COSOMG
      KHI = DASIN(SINKHI)
      Q = DSIN(OMEGA)
      R = COSOMG * DCOS(KHI)
      COSFI = (Q*T1+R*T2)/(T1*T1+T2*T2)
      SINFI = (R*T1-Q*T2)/(T1*T1+T2*T2)
      IF (SINFI .GT. -.7) GOTO 4
          FI = 0.0 - DACOS (COSFI)
          GOTO 999
4     IF (SINFI .LT. .7) GOTO 5
          FI = DACOS (COSFI)
          GOTO 999
5     IF (COSFI .GT. 0) GOTO 6
          FI = PI - DASIN (SINFI)
          GOTO 999
6     FI = DASIN(SINFI)
C
999   FI = FI*180./PI
      KHI = KHI*180./PI
      TWOT = TWOT*180./PI
C
998   OMEGA = OMEGA*180./PI
      RETURN
      END
